%% CLEAR ALL 
clear all;
close all;
clc;
rng(1, 'twister');
tic;

%% ADD PATH 
% Determine where your m-file's folder is.
folder = fileparts(which("Emp_App.m")); 
% Add that folder plus all subfolders to the path.
addpath(genpath(folder)); 
% Saving path
file_dir = what('Results'); 
spath=file_dir.path;
% Load Data
Z = xlsread('Data/Data.xlsx','raw data','B2:H203'); %load excel values
[Y,X,X_full,Q,Q_new]=transform(Z);


% Settings
hor=24;                            % Forecast horizons
jhor=[0:hor];                      % Time line
p=4;                               % Lag length
n_boot=500;                        % Number of bootstraps
weight=0.5;                        % Weighting for VAR IRFs
nP=p;
nH=hor+1;
BLP_settings;                      % Bayesian Local Projection settings

% Data preparation
qt=[];                             % Other controls
yt=Q_new(:,1:end);                 % Set of interested variables
xt=yt(:,end);                      % Federal funds rate
N_yt=size(yt,2);                   % Number of yt variables
RGDP_pos=1;  DGDP_pos=2;           % Positions of variables
Cons_pos=3;  GPDInv_pos=4;  
EmpHours_pos=5;	RealCompHour_pos=6;  
FedFunds_pos=7;
transformcode=Z(2,:);              % Transformation code for all variables
transformcode_varlp=Z(2,1:N_yt);   % Transformation code for interested variables

% Estimate IRFs by VAR
Y_var=Q(3:end,1:N_yt);
[Y_var_lagsnon1,Y_var_lags,Y_var_ad]=lag_variable(Y_var,p);
[T_var,N_var]=size(Y_var_ad);
[M_var]=size(Y_var_lags,2);

[IRF_VAR_full]=VAR_model(Y_var,hor,p,1,0);
IRF_VAR_org=IRF_VAR_full(end-N_var+1:end,:); 
IRF_VAR_result=IRF_VAR_org./IRF_VAR_org(end,1);
[IRF_VAR_convert]=IRF_convert(IRF_VAR_result,N_var,hor,transformcode_varlp);

% Estimate IRFs by BLP
IRF_BLP_full=IRF_BLP_estimation(Y_var,p,hor+1,Y_var_lagsnon1,Y_var_lags,Y_var_ad,hyperPars,hyperPriorsOptions,0);
IRF_BLP_org=IRF_BLP_full(end-N_var+1:end,:); 
IRF_BLP_result=IRF_BLP_org./IRF_BLP_org(end,1);
[IRF_BLP_convert]=IRF_convert(IRF_BLP_result,N_var,hor,transformcode_varlp);

% Estimate IRFs by LP
[IRF_LP_full]=LP_Jorda_model(Y_var,hor,p,1,0);
IRF_LP_org=IRF_LP_full(end-N_var+1:end,:); 
IRF_LP_result=IRF_LP_org./IRF_LP_org(end,1);
[IRF_LP_convert]=IRF_convert(IRF_LP_result,N_var,hor,transformcode_varlp);

% Sequence of block bootstrap
[~,~,~,~,IRF_VAR_boot_full,IRF_LP_boot_full]=VARLP_seqblockbootstrap(Y_var,p,hor+1,n_boot,1,1,0,0,0);
IRF_VAR_boot_org=IRF_VAR_boot_full(end-N_var+1:end,:,:); 
IRF_LP_boot_org=IRF_LP_boot_full(end-N_var+1:end,:,:); 

for i=1:n_boot
IRF_VAR_boot_convert(:,:,i)=IRF_convert(IRF_VAR_boot_org(:,:,i),N_var,hor,transformcode_varlp);
IRF_LP_boot_convert(:,:,i)=IRF_convert(IRF_LP_boot_org(:,:,i),N_var,hor,transformcode_varlp);
end

for i_n=1:N_var
for j=1:hor+1
covariance_varlp_seqblock(:,:,j,i_n)=cov(reshape(IRF_VAR_boot_convert(i_n,j,:),[n_boot,1,1]),reshape(IRF_LP_boot_convert(i_n,j,:),[n_boot,1,1]));
correlation_varlp_seqblock(:,:,j,i_n)=corrcoef(reshape(IRF_VAR_boot_convert(i_n,j,:),[n_boot,1,1]),reshape(IRF_LP_boot_convert(i_n,j,:),[n_boot,1,1]));
end
end

% Residual based bootstrap
[~,~,~,~,IRF_VAR_bootrep,IRF_LP_bootrep]=VAR_bootstrap_replacement(Y_var,hor,p,1,n_boot,1,1);
IRF_VAR_bootrep_org=IRF_VAR_bootrep(end-N_var+1:end,:,:); 
IRF_LP_bootrep_org=IRF_LP_bootrep(end-N_var+1:end,:,:); 

for i=1:n_boot
IRF_VAR_bootrep_convert(:,:,i)=IRF_convert(IRF_VAR_bootrep_org(:,:,i),N_var,hor,transformcode_varlp);
IRF_LP_bootrep_convert(:,:,i)=IRF_convert(IRF_LP_bootrep_org(:,:,i),N_var,hor,transformcode_varlp);
end

for i_n=1:N_var
for j=1:hor+1
covariance_varlp_bootrep(:,:,j,i_n)=cov(reshape(IRF_VAR_bootrep_convert(i_n,j,:),[n_boot,1,1]),reshape(IRF_LP_bootrep_convert(i_n,j,:),[n_boot,1,1]));
correlation_varlp_bootrep(:,:,j,i_n)=corrcoef(reshape(IRF_VAR_bootrep_convert(i_n,j,:),[n_boot,1,1]),reshape(IRF_LP_bootrep_convert(i_n,j,:),[n_boot,1,1]));
end
end

% Estimate IRFs by Pooling
IRF_Pool_convert=weight*IRF_VAR_convert+(1-weight)*IRF_LP_convert;


% Construct confidence intervals
scale_VAR=1/ IRF_VAR_org(end,1); 
scale_LP =1/ IRF_LP_org(end,1);
zvalue=1.96;

for i_n=1:N_var
for j=1:hor+1
IRF_VAR_variance_seqblock=covariance_varlp_seqblock(1,1,j,i_n);
IRF_LP_variance_seqblock=covariance_varlp_seqblock(2,2,j,i_n);
IRF_VARLP_covariance_seqblock=covariance_varlp_seqblock(2,1,j,i_n);

IRF_VAR_variance_bootrep=covariance_varlp_bootrep(1,1,j,i_n);
IRF_LP_variance_bootrep=covariance_varlp_bootrep(2,2,j,i_n);
IRF_VARLP_covariance_bootrep=covariance_varlp_bootrep(2,1,j,i_n);


% Variance of Pooled IRFs
IRF_Pool_sd_seqblock(i_n,j)=sqrt(weight^2*scale_VAR^2*IRF_VAR_variance_seqblock+(1-weight)^2*scale_LP^2*IRF_LP_variance_seqblock...
                              +2*weight*(1-weight)*scale_VAR*scale_LP*IRF_VARLP_covariance_seqblock);

IRF_Pool_sd_bootrep(i_n,j)=sqrt(weight^2*scale_VAR^2*IRF_VAR_variance_bootrep+(1-weight)^2*scale_LP^2*IRF_LP_variance_bootrep...
                              +2*weight*(1-weight)*scale_VAR*scale_LP*IRF_VARLP_covariance_bootrep);
end
end

% Confidence intervals
IRF_Pool_UB_seqblock=IRF_Pool_convert+zvalue*IRF_Pool_sd_seqblock;
IRF_Pool_LB_seqblock=IRF_Pool_convert-zvalue*IRF_Pool_sd_seqblock;

IRF_Pool_UB_bootrep=IRF_Pool_convert+zvalue*IRF_Pool_sd_bootrep;
IRF_Pool_LB_bootrep=IRF_Pool_convert-zvalue*IRF_Pool_sd_bootrep;

% Save results
Run_Emp_App_save_results
